home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 293_02 / grad1.c < prev    next >
C/C++ Source or Header  |  1989-08-23  |  33KB  |  1,012 lines
  1. /******************** GRAD ******************************
  2.   grad.c
  3.   Daniel Geist - Mallinckrodt Institue of Radiology 1987.
  4.    This program reads a set of ct scans for one
  5.    patient. It outputs two sets of files containing predefined 
  6.    views of bone surfaces. One set contains distance values and
  7.    the other gradient values.
  8.  
  9.    Danny Geist & Mike Vannier
  10.  
  11.    Contact for further information:
  12.     Michael W. Vannier
  13.     Mallinckrodt Institute of Radiology
  14.     Washington University School of Medicine
  15.     510 South Kingshighway Blvd.
  16.     St. Louis, Mo. 63110 USA
  17.  
  18.    Ake Wallin has made the following changes - 1988
  19.  
  20.    -Command line input
  21.    -Threshold of gradient for differencing
  22.    -Views to choose: BO (bottom), TO (top), RL (right lateral),
  23.     LL (left lateral), RE (rear), FR (front), NO (none, just leave distance)
  24.    -Clipping planes: RL (right lateral), LL (lateral left), RE (rear),
  25.     FR (front)
  26.    -Scan of data in "correct" order, assuming:
  27.     object:
  28.     low z - feet, high z - head
  29.     low x - right, high x - left
  30.     low y - back, high y - front
  31.     image:
  32.     low x - left, high x - right
  33.     low y - down, high y - up
  34.             front
  35.     view BO : right      left
  36.             back
  37.             front
  38.     view TO : left       right
  39.             back
  40.             head
  41.     view RL : back      front
  42.             feet
  43.             head
  44.     view LL : front     back
  45.             feet
  46.             head
  47.     view RE : left      right
  48.             feet
  49.             head
  50.     view FR : right     left
  51.             feet
  52. *********************************************************/
  53. #include <stdio.h>
  54. #include <math.h>
  55. #include <ctype.h>
  56.  
  57. /* some global variables*/
  58. int FIRSTSLICE,LASTSLICE,THRESHOLD,NLINES;
  59. float ZOOM;
  60. char DR;
  61. int  buffer[2][6][256];          /* input buffer */
  62. float fxbuf[5][256],fybuf[5][256]; /* X,Y Views floating buffers */
  63.  
  64. float huge fzbuf[256][256];        /* Z view floating buffer */
  65. char fnamein[50]="ctbild.000";        /* input file name  */
  66. float GRAD_THRESHOLD = 10.0;
  67. unsigned int   views = 0xFFFF;
  68. int   clipx[2], clipy[2];      /* clipping planes */
  69. int   dispmode;                /* distance and/or gradient */
  70. int   image_or[2];            /* image oriention in x and y */
  71. int   object_or[3];           /* image oriention in x, y and z */
  72. int   header_blocks;          /* number of header blocks in CT's */
  73. char *usestr1 = "Usage: grad [file] [-z] [-f] [-l] [-t] [-d] [-n(g|d)] [-g]";
  74. char *usestr2 = "            [-h] [-v(bo|to|rl|ll|re|fr|no)] [-c(rl|ll|re|fr)]";
  75. char *usestr3 = "            [-i(x(r|l)|y(u|d))] [-o(x(r|l)|y(f|b)|z(h|f))]";
  76.  
  77. succ(i)   /* (i+1) modolus 3 */
  78. int i;
  79. {
  80.     return(i==2?0:i+1);
  81. }
  82.  
  83. prev(i)   /* (i-1) modolus 3 */
  84. int i;
  85. {
  86.     return(i==0?2:i-1);
  87. }
  88.  
  89. /* Set input file - add slice number to extension */
  90. setfilename(filenum)
  91. int filenum;
  92. {
  93.     int i;
  94.  
  95.     for (i = strlen(fnamein)-1; i > 0; i--)
  96.       if (fnamein[i] == '.') break;
  97.     if (i == 0) {
  98.       printf("CT file names inconsistend!\n");
  99.       exit(1);
  100.     }
  101.     fnamein[++i]=filenum/100+'0';
  102.     fnamein[++i]=(filenum%100)/10+'0';
  103.     fnamein[++i]='0'+filenum%10;
  104. }
  105.  
  106. /* interpolate from bottom line to top line n-1 lines  */
  107. interpolate(line,bot,top,n)
  108. int line,bot,top,n;
  109. {
  110.     int next,i,j,x,inc;
  111.     inc=top>bot?1:(-1);    /* interpolate backward or forwards ? */
  112.     next=bot+inc;
  113.     for(i=1,j=n-1;i<n;i++,j--){   /* do for each next line of interpolation */
  114.         for(x=0;x<256;x++) buffer[line][next][x]=
  115.             (buffer[line][bot][x]*j+buffer[line][top][x]*i)/n;
  116.         next+=inc;
  117.     }
  118. }
  119.  
  120. /* midpoint - fraction part of threshold transition distance */ 
  121. float midpoint(b,a)
  122. float b,a;
  123. {
  124.     return( (THRESHOLD-a) / (b-a) );
  125. }
  126.  
  127. /* get floating point distance values  */
  128. getdistances(xstart,xstop,xdir,
  129.              ystart,ystop,ydir,
  130.              start_slice,end_slice,zdir,
  131.              xbufxstart, xbufxdir, xbufystart, xbufydir,
  132.              ybufxstart, ybufxdir, ybufystart, ybufydir,
  133.              zbufxstart, zbufxdir, zbufystart, zbufydir,
  134.              pass)
  135. int xstart,xstop,xdir,ystart,ystop,ydir,start_slice,end_slice,zdir,pass;
  136. int  xbufxstart, xbufxdir, xbufystart, xbufydir;
  137. int  ybufxstart, ybufxdir, ybufystart, ybufydir;
  138. int  zbufxstart, zbufxdir, zbufystart, zbufydir;
  139. /* also used are global variables: fxbuf, fybuf, fzbuf, buffer,
  140.    ZOOM, DR, image_or, object_or. */
  141. {
  142.     int z,x,y,i,j,start,stop,inc,line,rzoom,inter;
  143.     float remain;
  144.     FILE *fxfloat,*fyfloat,*fzfloat,*fn[2];
  145.     char filename[13];
  146.     int  filegap;
  147.     int  xbfx, ybfx, zbfx, zbfy;
  148.  
  149.     NLINES=0;         /* number of output lines in X,Y view directions */
  150.     remain=0;         /* remainder of interpolation after roundoff */
  151.     rzoom=ZOOM+0.5;   /* rounded zoom factor */
  152.     /* X view floating output file */
  153.     sprintf(filename,"%c:xdis%d.dat",DR,pass);
  154.     if ((fxfloat=fopen(filename,"wb")) == NULL) {
  155.       printf("Error creating %s!\n", filename);
  156.       exit(1);
  157.     }
  158.     /* Y view floating output file */
  159.     sprintf(filename,"%c:ydis%d.dat",DR,pass);
  160.     if ((fyfloat=fopen(filename,"wb")) == NULL) {
  161.       printf("Error creating %s!\n", filename);
  162.       exit(1);
  163.     }
  164.     for(i=0;i<256;i++)
  165.         for(j=0;j<256;j++)fzbuf[i][j]=256; /* clear Z view buffer */
  166.     
  167.     for(z=start_slice; z!=end_slice; z += zdir){  /* For each Slice */
  168.       /* open next two slice files */
  169.       setfilename(z);
  170.       if ((fn[0]=fopen(fnamein,"rb")) == NULL) continue;
  171.     for (filegap = zdir; filegap + z != end_slice+zdir; filegap+=zdir) {
  172.           setfilename(z+filegap);
  173.           if ((fn[1]=fopen(fnamein,"rb")) != NULL) break;
  174.       }
  175.       if (fn[1] == NULL) continue;
  176.  
  177.       inter=rzoom; /* interpolation factor assumed rounded zoom */
  178.       /* correct interpolation factor according to floating remainder */
  179.       remain+=rzoom-ZOOM;
  180.       if(remain>=1){
  181.         inter-=1;
  182.         remain-=1;
  183.       }
  184.       else if(remain<=(-1)){
  185.         inter+=1;
  186.         remain+=1;
  187.       }
  188.  
  189.       line=0;               /* current input buffer line */
  190.       for(j=0;j<inter;j++)  /* clear X,Y floating buffers */
  191.         for(i=0;i<256;i++) fxbuf[j][i]=fybuf[j][i]=256;
  192.       /* set index for zbuf */
  193.       zbfy = zbufystart;
  194.       /* set index for xbuf */
  195.       xbfx = xbufxstart;
  196.       for(y=ystart; y!=ystop; y+=ydir){              /* For each line */
  197.     /* skip to line*/
  198.     for(i=0;i<2;i++)fseek(fn[i],(long)512*(y+header_blocks),SEEK_SET);
  199.         fread(buffer[line][0],1,512,fn[0]);   
  200.         fread(buffer[line][inter],1,512,fn[1]);
  201.         interpolate(line,0,inter,inter); /* interpolate in_between */
  202.     for(i=0;i<inter;i++){   /* For each interpolation line */
  203.           /* set index for zbuf */
  204.       zbfx = zbufxstart;
  205.           /* set index for ybuf */
  206.           ybfx = ybufxstart;
  207.       for(x=xstart; x!=xstop; x+=xdir) { /* For each Voxel value */
  208.             /* find threshold transition*/
  209.         if (buffer[line][i+1][x] >= THRESHOLD) {
  210.               /* if first transition in X direction get floating
  211.                  distance */
  212.               if(fxbuf[i][xbfx]==256.0) fxbuf[i][xbfx]=(x==xstart)?0:
  213.                                     (x-xstart)*xdir-1+
  214.                      midpoint((float)buffer[line][i+1][x],
  215.                             (float)buffer[line][i+1][x-xdir]);
  216.  
  217.               /* if first transition in Y direction get floating
  218.                    distance */
  219.               if(fybuf[i][ybfx]==256.0) fybuf[i][ybfx]=(y==ystart)?0:
  220.                     (y-ystart)*ydir-1+
  221.                      midpoint((float)buffer[line][i+1][x],
  222.                               (float)buffer[1-line][i+1][x]);
  223.  
  224.               /* if first transition in Z direction get floating
  225.                    distance */
  226.               if(fzbuf[zbfy][zbfx]==256.0) fzbuf[zbfy][zbfx]=
  227.                    (i==0) && (buffer[line][i][x]>=THRESHOLD) ?
  228.                     NLINES : NLINES+i+
  229.                     midpoint((float)buffer[line][i+1][x],
  230.                              (float)buffer[line][i][x]);
  231.             }
  232.             /* change index of ybuf and zbuf */
  233.             ybfx += ybufxdir;
  234.         zbfx += zbufxdir;
  235.         if ((ybfx < 0 || ybfx >